Model 3: SEM CM without respecifications

Results

##   Soil_property                    ME              RMSE               SS
## 1         CEC.A   0.00775128089893602  5.68793208043344 5047.00113085333
## 2         CEC.B  -0.00398182488814064  6.26155391338747 6116.30095600024
## 3         CEC.C   0.00229007311875548  6.77752944249498 7165.84523364626
## 4          OC.A  0.000480842199394078 0.629623927110415 61.8425011760308
## 5          OC.B  0.000916586572794572 0.316970003100663  15.673317327039
## 6          OC.C  0.000757381666383179 0.180989934293459 5.11014758522591
## 7        clay.A   0.00387041702804282  8.34581735151287 10865.8160933109
## 8        clay.B   -0.0125981787909463  9.39142974125339 13759.0366032442
## 9        clay.C -0.000429563387851167  10.0827907785092 15859.3765017776
##   mean_theta median_th         R2
## 1  0.9450710 0.1988536 0.05396970
## 2  0.7638961 0.2693472 0.15050286
## 3  0.7684111 0.3220226 0.16161086
## 4  1.0386352 0.3081390 0.06873298
## 5  0.9604539 0.2237897 0.02997153
## 6  1.4402485 0.1209572 0.04206266
## 7  0.8671333 0.1784460 0.11736397
## 8  0.8878465 0.3450671 0.21954429
## 9  1.1430615 0.2243374 0.15005788
##   Soil_property            ME     RMSE        r2
## 1          CECo  0.0020198430 6.258182 0.2025248
## 3           OCo  0.0007182701 0.420180 0.6836467
## 5         clayo -0.0030524417 9.300794 0.2659860

Spatial autocorrelation analysis

### Autocorrelation in SEM residuals ####
backup <- Res
Res[2:10] <- Res[,2:10]-Res[,11:19]
Res <- Res[,1:10]
names(Res)[1] <- "id.p"

R <- merge(Res, unique(d[,c(1,20,21)]), by.x = "id.p", by.y = "idp", all.x=T)
library(sp)
library(gstat)
# R as geo
coordinates(R) <- ~X+Y
#define crs
wgs84 <- CRS("+init=epsg:4326")
#UTM14N <- CRS("+init=epsg:32614")
modis <- CRS("+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs")
NAD83.KS.N <- CRS("+init=epsg:2796")

proj4string(R) <- NAD83.KS.N
R <- spTransform(R, wgs84)
 mapview::viewExtent(R)
library(mapview)
## Loading required package: leaflet
mapview(R, burst = TRUE)
###############################
# Variograms
# independent of device size
g <- list()
vg <- list()
vgm <- list()
# CEC
for(i in 2:4){
  g[[i-1]] <-  gstat(id = c(names(R@data)[i]), 
                     formula = formula(paste0(names(R@data)[i],"~1")),
            data = R)
  vg[[i-1]] <-  variogram(g[[i-1]])
  # vg = variogram(g, width = 20000, cutoff = 600000)
  # vg = variogram(g, boundaries = c(1E4,2E4,4E4,7E4,1E5,2E5,4E5,7E5,1E6))
  # print(plot(vg, plot.numbers = TRUE, main = names(R@data)[i]))
  
  # # choose initial variogram model and plot:
  vgm[[i-1]] <- vgm(nugget = var(R@data[,i]),
                    psill=var(R@data[,i])*2,
                    range=100000,
                    model = "Sph")
  #plot(vg, vgm, main = names(R@data)[i])
  # 
  # # fit variogram model:
  vgm[[i-1]] <-  fit.variogram(vg[[i-1]], vgm[[i-1]])
  print(names(R@data)[i])
  print(vgm[[i-1]])
  print(attr(vgm[[i-1]], "SSErr"))
}
## Warning in fit.variogram(vg[[i - 1]], vgm[[i - 1]]): No convergence after
## 200 iterations: try different initial values?
## [1] "CEC.A"
##   model    psill    range
## 1   Nug 32.56124      0.0
## 2   Sph 65.12263 100000.2
## [1] 487.1739
## Warning in fit.variogram(vg[[i - 1]], vgm[[i - 1]]): No convergence after
## 200 iterations: try different initial values?
## [1] "CEC.B"
##   model    psill   range
## 1   Nug 39.45999     0.0
## 2   Sph 78.91982 99999.8
## [1] 1681.899
## Warning in fit.variogram(vg[[i - 1]], vgm[[i - 1]]): No convergence after
## 200 iterations: try different initial values?
## [1] "CEC.C"
##   model    psill    range
## 1   Nug 46.23125     0.00
## 2   Sph 92.46235 99999.83
## [1] 4387.061
# OC
for(i in 5:7){
  g[[i-1]] <-  gstat(id = c(names(R@data)[i]), 
                     formula = formula(paste0(names(R@data)[i],"~1")),
            data = R)
  vg[[i-1]] <-  variogram(g[[i-1]])
  # vg = variogram(g, width = 20000, cutoff = 600000)
  # vg = variogram(g, boundaries = c(1E4,2E4,4E4,7E4,1E5,2E5,4E5,7E5,1E6))
  # print(plot(vg, plot.numbers = TRUE, main = names(R@data)[i]))
  
  # # choose initial variogram model and plot:
  vgm[[i-1]] <- vgm(nugget = var(R@data[,i]),
                    psill= var(R@data[,i])*2,
                    range=100000,
                    model = "Sph")
  #plot(vg, vgm, main = names(R@data)[i])
  # 
  # # fit variogram model:
  vgm[[i-1]] <-  fit.variogram(vg[[i-1]], vgm[[i-1]])#, fit.method = 4)
  print(names(R@data)[i])
  print(vgm[[i-1]])
  print(attr(vgm[[i-1]], "SSErr"))
}
## Warning in fit.variogram(vg[[i - 1]], vgm[[i - 1]]): No convergence after
## 200 iterations: try different initial values?
## [1] "OC.A"
##   model     psill   range
## 1   Nug 0.3989836     0.0
## 2   Sph 0.7979665 99999.9
## [1] 0.1327332
## Warning in fit.variogram(vg[[i - 1]], vgm[[i - 1]]): No convergence after
## 200 iterations: try different initial values?
## [1] "OC.B"
##   model     psill    range
## 1   Nug 0.1011173      0.0
## 2   Sph 0.2022348 100000.1
## [1] 0.008639608
## Warning in fit.variogram(vg[[i - 1]], vgm[[i - 1]]): No convergence after
## 200 iterations: try different initial values?
## [1] "OC.C"
##   model      psill    range
## 1   Nug 0.03296812      0.0
## 2   Sph 0.06593627 100000.1
## [1] 0.0009133788
# Clay
for(i in 8:10){
  g[[i-1]] <-  gstat(id = c(names(R@data)[i]), 
                     formula = formula(paste0(names(R@data)[i],"~1")),
                     data = R)
  vg[[i-1]] <-  variogram(g[[i-1]])
  # vg = variogram(g, width = 20000, cutoff = 600000)
  # vg = variogram(g, boundaries = c(1E4,2E4,4E4,7E4,1E5,2E5,4E5,7E5,1E6))
  # print(plot(vg, plot.numbers = TRUE, main = names(R@data)[i]))
  
  # # choose initial variogram model and plot:
  vgm[[i-1]] <- vgm(nugget = var(R@data[,i]),
                    psill=var(R@data[,i])*2,
                    range=5E4,
                    model = "Sph")
  #plot(vg, vgm, main = names(R@data)[i])
  # 
  # # fit variogram model:
  vgm[[i-1]] <-  fit.variogram(vg[[i-1]], vgm[[i-1]], fit.method = 2)
  print(names(R@data)[i])
  print(vgm[[i-1]])
  print(attr(vgm[[i-1]], "SSErr"))
}
## Warning in fit.variogram(vg[[i - 1]], vgm[[i - 1]], fit.method = 2): No
## convergence after 200 iterations: try different initial values?
## [1] "clay.A"
##   model     psill   range
## 1   Nug  70.10202     0.0
## 2   Sph 140.20911 50001.8
## [1] 336.8401
## Warning in fit.variogram(vg[[i - 1]], vgm[[i - 1]], fit.method = 2): No
## convergence after 200 iterations: try different initial values?
## [1] "clay.B"
##   model     psill    range
## 1   Nug  88.76782     0.00
## 2   Sph 177.53011 49998.44
## [1] 299.7812
## Warning in fit.variogram(vg[[i - 1]], vgm[[i - 1]], fit.method = 2): No
## convergence after 200 iterations: try different initial values?
## [1] "clay.C"
##   model    psill    range
## 1   Nug 102.3186     0.00
## 2   Sph 204.6411 50000.97
## [1] 621.602
# Three graphs (soil properties)
par(mfrow = c(1, 3), pty = "s")
### CEC
## A
plot(variogramLine(vgm[[1]], maxdist=2E5), 
     type="l", lwd=2,col="#AA0000", main="CEC",
     xlab = "Distance", ylab = "Semivariance", cex.lab = 1.3, ylim=c(0,50))
points(gamma ~ dist, vg[[1]], col="#770000")
#legend(x= "topleft",legend = "SSErr", bty = "n")
#legend(x= 12,legend = round(attr(vgm[[1]], "SSErr"),3), text.col ="#AA0000", 
#        bty = "n")
## B
lines(variogramLine(vgm[[2]], maxdist=2E5), lwd=2, col="#00AA00")
points(gamma ~ dist, vg[[2]], col="#007700")
## C
lines(variogramLine(vgm[[3]], maxdist=2E5), lwd=2, col="#0000AA")
points(gamma ~ dist, vg[[3]], col="#000077")

### OC
## A
plot(variogramLine(vgm[[4]], maxdist=2E5), type="l", lwd=2,col="#AA0000", 
     main="OC",
     xlab = "Distance", ylab = "Semivariance", cex.lab = 1.3, ylim=c(0,0.5)) 
points(gamma ~ dist, vg[[4]], col="#770000")
## B
lines(variogramLine(vgm[[5]], maxdist=2E5), lwd=2, col="#00AA00")
points(gamma ~ dist, vg[[5]], col="#007700")
## C
lines(variogramLine(vgm[[6]], maxdist=2E5), lwd=2, col="#0000AA")
points(gamma ~ dist, vg[[6]], col="#000077")

### Clay
## A
plot(variogramLine(vgm[[7]], maxdist=2E5), type="l", lwd=2,col="#AA0000", 
     main="Clay",
     xlab = "Distance", ylab = "Semivariance", cex.lab = 1.3, ylim=c(0,140)) 
points(gamma ~ dist, vg[[7]], col="#770000")
## B
lines(variogramLine(vgm[[8]], maxdist=2E5), lwd=2, col="#00AA00")
points(gamma ~ dist, vg[[8]], col="#007700")
## C
lines(variogramLine(vgm[[9]], maxdist=2E5), lwd=2, col="#0000AA")
points(gamma ~ dist, vg[[9]], col="#000077")
#############################################################################
par(mfrow = c(1, 1))

title("Semivariograms of residuals",
  sub="red: A horizon, green: B horizon, blue: C horizon")

# vgm[[4]] <- vgm(nugget = 0.1,
#                   psill=0.1,
#                   range=5E4,
#                   model = "Wav")
# plot(vg[[4]], vgm[[4]],ylim=c(0,0.5))
# vgm[[4]] <-  fit.variogram(vg[[4]], vgm[[4]], fit.method = 2)
# plot(vg[[4]], vgm[[4]],ylim=c(0,0.5))
# attr(vgm[[4]], "SSErr")
# vgm[[4]]